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Целью данной работы является рассмотрение пространствен- 
но-трехмерной модели Баклея — Леверетта фильтрации 
двухфазной несжимаемой жидкости «нефть — вода» в при- 
ближении Буссинеска и ее численная реализация усовершен- 
ствованным вариантом модифицированного попеременно- 
треугольного метода (МПТМ), учитывающим специфику 
функций точечных источников, имеющих дельтаобразный 
характер. Для предложенного усовершенствованного варианта 
МПТМ выполнено построение набора итерационных пара- 
метров и получены необходимые оценки спектральных пара- 
метров метода. Оценка скорости сходимости итерационного 
процесса асимптотически эквивалентна оценке для базового 
итерационного МПМ, предложенного ранее 
А. А. Самарским и Е. С. Николаевым. Построенный вариант 
итерационного метода был применен для решения модельной 
задачи фильтрации — вытеснения нефти методом заводнения 
в существенно неоднородном пласте, для которого значения 
коэффициента проницаемости менялись более чем в 104 раза. 
Что касается вычислительной эффективности такого подхода, 
то число итераций сократилось в 10—70 раз по сравнению с 
традиционно используемыми итерационными методами — 
Зейделя и верхней релаксации с шахматной упорядоченно- 
стью узлов. 


Ключевые слова: модели фильтрации Баклея — Леверетта, 
двухфазная несжимаемая жидкость, разностные схемы, итера- 
ционный попеременно-треугольный метод, итерационные 
параметры. 
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Введение. В современных условиях для решения значимых задач проектирования разработок нефтяных ме- 
сторождений (РНМ) и совершенствования технологии добычи углеводородного сырья используется методология ма- 
тематического моделирования и комплексы программ, реализующих современные вычислительные методы [1, 2]. 
Важный класс задач, возникающих при проектировании РНМ, — задачи фильтрации многофазных жидкостей в пори- 
стых средах, сформулированные при тех или иных упрощающих предположениях. Данный класс задач после дискре- 
тизации приводит к сеточным уравнениям, которые могут иметь размерности 10—10” и должны решаться многократ- 
но (многие десятки, сотни тысяч раз) для временных промежутков, составляющих десятки месяцев. В практике вы- 
числительных методов решения задач данного класса широко применяется /МРЕ5-метод (ПМрйси оп Рте5биге апа Ех- 
РИсй оп Зашганоп) — неявный по давлению и явный по насыщенности [3]. Поэтому основной объем вычислительной 
работы приходится на решение уравнений эллиптического типа в относительно небольшом числе узлов сетки с функ- 
циями источников, отличными от нуля и имеющими характер дельта-функций (функций точечных источников) для 
определения давления. Следует учитывать данную специфику сеточного оператора задачи. В противном случае, а 
также при недостаточной обусловленности системы разностных уравнений, вызванной высокой размерностью и су- 
щественной неоднородностью характеристик пластовой системы, снижается качество вычислений. В первую очередь 
следует отметить, что замедляются итерационные процессы, а при использовании техники чебышевского ускорения 
наблюдается потеря сходимости метода, обусловленная ошибкой в определении спектральных характеристик опера- 
торов. В данной работе указанная проблема снимается за счет построения оператора-предобусловливателя, содержа- 
щего член, обусловленный наличием источников и корректных спектральных двухсторонних оценок для оператора- 
предобусловливателя. 

1. Непрерывная пространственно-трехмерная модель двухфазной фильтрации. Нефтеводоносный пласт 
может существенным образом менять глубину залегания и (или) иметь существенную толщину и неоднородность по 
вертикали. В таких случаях требуется использование пространственно-трехмерной модели, учитывающей влияние 
силы тяжести, толщины пласта и неоднородности его параметров по вертикали. Будем использовать Декартову систе- 
му координат — ось Ох направлена на север, ось Оу — на восток, Ох — вертикально вниз. Процесс фильтрации 


рассматривается в пространственно-трехмерной области С — цилиндрической области с боковой поверхностью 


> (х, У), (х, У) ЕГ, Г=В\ р, ограниченной снизу поверхностью подошвы пласта, задаваемой функцией поверхности 


&=Нь(х, у), (х,у)ЕБ,Ое В’, асверху — поверхностью кровли пласта = Н, (х,у), (х, у) ЕС. 

В отношении граничных поверхностей области С предполагается их достаточная гладкость, необходимая для 
корректной постановки задачи. Заметим, что в этом случае граничная поверхность Г области С состоит из объедине- 
ния достаточно гладких поверхностей Г=Н,(]Н,ЦУ, @в=Сс0Г. 


Для проницаемой кровли пласта уравнения, описывающие процесс фильтрации в переменных «давление — 
водонасыщенность», в случае несжимаемости сред имеют вид: 


9 | ЛСК, об Ор ы 9 ($), ОА Ор 





9х И и, 0х) ду и и бу 
в по [226% [2+ вр) а. (0 
05\ м 0% > 0 


т — д (5) К, Ор + д Ь (5)К, Ор 
9 Ох) 2 9х) 9) 12 9 


д к, (д 
ЕВА 26 [+ вр, + 42. (2) 
05\ № \0 


В уравнениях (1), (2): 5 = 5(х, у, 5,1) — водонасыщенность; } (5), }›(5) — относительные фазовые проницае- 





МОСТИ ДЛЯ нефти и воды соответственно; и — пористость пласта; д, #› — динамические вязкости нефти и воды С0- 
ответственно; р, р› —— плотности нефти и воды соответственно, а ру, р› — их фазовые давления соответственно; 
К ‚ К, = коэффициенты абсолютной проницаемости соответственно в горизонтальном и вертикальном направлениях. 
Здесь мы отметим следующее: как правило, К, значительно больше К, и К, существенно изменяется в зависимости 


от координаты с. 
Будем считать, что на скважинах задаются либо дебиты, либо забойные давления. В качестве функциональных 


зависимостей для задания Л (5), = 1, 2 будем использовать полиномы второго порядка относительно переменной 5. 
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Граничные условия будем рассматривать в связи с тем, является ли граница Г области С проницаемой. На 


ь д В 

непроницаемои границе поток по нормали должен быть равен нулю, что приводит к условию — =0. На проницаемои 
п 

границе, которой также является контур, ограничивающий скважину, будем рассматривать граничные условия 1-го и 


2-го рода. Пусть и, — суммарный поток многофазной жидкости. Предположим, совместное движение фаз 
И, = И +, 


2$) 
где и," — потоки нефти и воды соответственно, удовлетворяющие условиям: у =-—К- ^^ ота4 р. 


1 

Граничные условия 2-го рода реализуются, когда задан поток, закачиваемый (отбираемый) на части цилин- 

дрической границы %/: 
и’, (х, у, 2,1) =М (х, у,7,1), (хузеЕХ, 0<Е5Цц, 
5(х, у, 7,1) = 5 (х, у, 2,1), (ху) ЕХ, 0<Е<Цц, 
а также для отбора, заданного на части цилиндрической границы Х.: 
и’, (х, у, 1) =У/(х, у,1), (ху) ЕХ,, О<Е5Ц, 
или давления 
и’, (х, УВ) =У(х, у,Ю), (х, у, 5) Е Хо, О<Е3ц, 
Р(х, у.) = Р(х,у,1), (ху ЕХ., О<Е5Щ. 

Разумеется, на части цилиндрической поверхности У. также могут быть заданы условия непроницаемости. 

2. Дискретизация пространственно-трехмерной модели фильтрации на основе метода неявного по дав- 
лению и явного по насыщенности. Ранее МПТМ для двумерных задач фильтрации с огрубленными спектральными 
оценками рассматривался в работе [4], а для уравнений параболического типа — в монографии [5]. МПТМ для случая 
несамосопряженного оператора сеточной задачи при ограниченном значении числа Пекле, предназначенный для ре- 
шения задач морской гидродинамики, был построен и исследован в работах [6, 7]. Далее рассматривается как наибо- 
лее распространенный алгоритм дискретизации задачи, базирующийся на [МРЕ$ — методе, неявном по давлению и 
явном по насыщенности. Также будем считать, что С(х, у, <) — область, в которой требуется найти решение задачи 


(1)—(2), — параллелепипед. В области С используем равномерную пространственную сетку @, : 





@, ={ ж=М,, у, = Т,, ак =#1,,0<1<М,,0<]<М,,0<<М,} 
и неравномерную временную сетку @,. 
Аналогично предыдущему определим пространственно-временную сетку (2) = [2 х@, В 


Разностные аппроксимации уравнений (1) и (2) во внутренних узлах сетки @, имеют вид для уравнения (2): 


1 @ 1 (1) 
У (А-ьл Е у (В НИ у + 
х 2 Е 2 
г чье вы: 80 (Ви - Ва) В + 
Е р 
У 2 У 2 
(В еб 8% гв он] В 44 
и2 ь.К-1 К лк и? К .К Бе 


я (25 (Арль + ЛЕ ков —0,5 (лк + ЛЕ )*, ЬрК+0,5 ) + 
а 





+ —. (0,5 (л и р-н —0,5 (ь ва Ля )к, Ь.К+0,5 ) — 


0 — если узел вне скважины; 


Воть Рурк В: | - еСли узел находится на нагнетательной скважине; 
8), "Е, ], 1,» 


(3) 


-В..,.Р..-Р. 2 
ВК оБрк Е }— если узел находится на эксплуатационной скважине. 


Вестник Донского государственного технического университета 2016, №1(84), 5-18 





Следует отметить, что значения функций относительной фазовой проницаемости в уравнении (3) и в уравне- 


нии (4) в зависимости от значений функции водонасыщенности вычисляются на предыдущем временном слое (метод 
явный по функции водонасыщенности): 

















п-+1 п 
511-55; , 1 
БУ ТЫ = + 
тН - 2 (5 Баре +9 РЕ, (ель в )+ 
п-1 ПХ 
1 = + 
+— 5 (5 Бык + 92 Фнльйннилк | (Вл И. )- 
— х 
2 (5 РрракКвы рык + 93 * льна (Р Би В.к)+ 
шв? 
1 + 
5 (54 АЙ + 64 Фаны, ны (В. -Руми)- 
мт 
— 
С (5% Лика +95 Локи, &)(в. Вл) + 
и, 
= ПР +95. Фльйи, ыы (В -Р кн) 
” < 
_ 822. 82 
ий, 222 (5; Ь: ‚К- я ‚Ьк— +55 Ь: „были )+ |’ ЙА 272 (5 Ь: РЕЧКУ + д5 Ь, ЕКы, к} = 
ГИ 2% 


0—если узел вне скважины; 


= Воть К-Рук - если узел находится на нагнетательной скважине; 
81,1, Г, Г], Л, 








(4) 
-Р.,5П., . 
1, К и к К)” ЬЛК - если узел находится на эксплуатационной скважине. 
Забойные давления на нагнетательной и эксплуатационной скважинах соответственно: 
Вр ЮБК: 
и (5) № (5) _ 2л 
Веру = АК кВ и Ей т . (5) 
ое Е: Ш— -15щ2-с 
та 
+ сос с с сх < 
1 ‚02, 03,04 ,05,665 9,02, 0,04 ,05 , 06 ы 
Символы определяются из условий: 
=(1-5, ), 6 =1,2,3,4,5,6 
хе 1-если В к -В, к 20, — 1-если В кВа: к 20, 
| = 5 = 
0 — в противном случае; 0 — в противном случае; 
5- 1- если В к -В, к 20, 5- 1- если В к -В, к 20, 
— 4 = 
0 — в противном случае; 0 — в противном случае; (6) 
зе 1-если В 1 -В, к 20, а 1- если В к-В кн 20, 
5 = 6 = 
0 — в противном случае; 0 — в противном случае. 


Приведем соотношения, определяющие коэффициенты разностных схем (3), (4), полученные интегрально- 
интерполяционным методом. 


--1 








Во В Дфах 
. 1 $ ыы 
НЫ | № х Ур (м (565 У») + У») 
, 1 
В©) т шр>4У 
Но п, у Ка (жь, У к) (ЛЬ, У, к) + Л У, 2%) 
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—1 


1 “= а 
а | Г 
р са вру (а (6,9) Ур в) 


3. Усовершенствование модифицированного попеременно-треугольного метода решения простран- 
ственно-трехмерной разностной задачи двухфазной фильтрации. Будем для простоты рассматривать граничные 
условия 1-го рода. Систему разностных уравнений (6) запишем в стандартном виде [8]: 


з 
-2 (= '›», +а(х)у(х)=Ф(х), хео, (8) 


х 
У(х)=и(х), хет, уеГ, (9) 
О<а о а(х)>0, Ухео. 
Здесь Г — граница прямоугольного параллелепипеда С = {0 <х, <, а=Ь2, 3} ‚0=®-+‘у — равномерная сетка, по- 
крывающая область С, у — множество граничных узлов сетки. 
Связь между коэффициентами уравнений (8) и (3) задается очевидным образом 


+1 _ _ р > — ра) 
а =а (м; +1, ;, 2х) = В 1 ра = (мрак) = В Г 
Н»/, = 


‚К 
59 
м_ — 5 Ц_ Ве, 
а> а, (©; +в, жк ) = ВО 1 ‚а =а> (хь› жж к ) = ВО 1 › 
1] +-К 1,]+-,К 
2 2 
+1 — — р(3) — _ р(3) 
аз = аз (жир ж жк +) = ВО -@ = 43 (мьж» жк} =В9 1) 
+ Би 


Хх Е @, Хх) Е @,, ЖЕ, 
@, ={х, =1й,,15<М,-Ц, а=Ъ2,3, Ф@=@ хо, ха. 


Предполагается, что имеет место связь между искомыми сеточными функциями уравнений (10) и (11) и се- 
точными функциями правых частей этих же уравнений: 


У(жиь м» жыь ) = И,» (мех) = Венло 





__ 82 
Ф(жи> хо хзк = ВР [5 ("ик о к, Вк-0,5 — 
й, 
8Р 
-0,5 (Л ильн + лк ре, лкноз) = т т (6,5(", вк Л )*, Бк-05 
2% 


п п 
-0,5(Х КО ры льнов | и 
Здесь значения относительных фазовых проницаемостей: /, ; а 1"; входят в соот- 
11, .к-Т› ЮУ 2ЬК-У 2 
ветствии с равенствами (10) в функции правых частей сеточного уравнения (8) и рассчитываются на основе значений 


сеточной функции водонасыщенности на предыдущем п-м временном слое. При этом 


В 7 р к если узел находится на нагнетательной скважине; 
> 


Р — если узел находится на эксплуатационной скважине; 


Е. 
81, К 91, 
0О-если узел не совпадает ни с одной из скважин. 
Будем использовать смещенные сетки 
+ _ Е : + _ —__# г. + _ + . 
© =[* =АА,1<1 < М}, ={х, =Ьй>,1<Ь <М, }, ол ={% =5/№,1< 5 <М,} 
и скалярные произведения 


(и) м (и,ъ) = У и(х)у(х) в , 


@ 


хе’ 
ыы = Х ибоьдь, (ии = Хх ибобдь, 
С № е@> хе" 


(оды = Х ибо, (или = Х и(убь, 


+ 
№ е@з № Е@; 
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(му = ХХ У и(х)ь(х)вЬь, 


ЖЕа»Х №4 ХЕ 
(о -х хх иены, (из, = ХХХ и(ауедыььь. 

дей х›е@» № е@3 мЕ® хе@, уе 

Представим систему разностных уравнений (8), (9) в виде операторного уравнения с однородными граничны- 
ми условиями, изменив соответствующим образом правую часть: 


о 
Ау=и, 


Ха 


о 3 о о 
45-3 [42 5. | +4У, хеЕ@, (10) 
а=1 < 


о 
У(х)=0, хеу. (11) 
Представим схему итерационного двухслойного модифицированного попеременно-треугольного метода в 


операторном виде: 


п-1 п 
(р+ев) р (р+ 8, + Ау" =, 
Тит (12) 





где операторы, входящие в уравнение (12), запишутся 


3 [а Е 
в [м . с у+ °} 
а=1 Г 











о 2, 6 (13) 
Е ан ад 1 
Ву=-х |, + у--9 | 
а=\ и 2, 6 (14) 


Видно, что операторы, введенные в соответствии с равенствами (13) и (14), являются сопряженными на мно- 


жестве сеточных функций, обращающихся в ноль на границе сетки в соответствии с равенством (10), т. е.: 
* 
К =-Е>. 
Очевидно, что 
А=А +К.. 


Получим необходимые оценки констант, входящих в операторные неравенства [8]: 
_ д 
50 < В +8,, ВО "< (В+К,). (15) 


Сеточная функция 4(х), определяющая элементы диагональной матрицы ДО, будет построена ниже. 
Рассмотрим соответствующие скалярные произведения: 


оо 13 а о ао | 93 а о ао | о 
ВЮ, у,у |= в. 2 2 | а _ 
| ру >] [18 Е ря ры он, > 697 


а а 











+1 2 о 
1 3|а о Чо |9 
= У | +—=“ 1]. 
а Г И, У», 21, . 6 т 
Будем использовать лемму: «Для любых чисел ри, Иа, Уи, Г. имеет место следующее неравенство: 
п 2 п тп 2 2 Я. 
>. [Рио +9ит, | 3 >. > вов ([ра|-+ Жа вена + м) (17) 
ед @= = Ха 


ГД &в > 0,58 = ри 20 = 9, > 0», 


Доказательство этой леммы основывается на =-неравенстве и здесь не приводится [7]. 
Положим в неравенстве (17) 
+1 




















а, . 1 р Ч йа 
==, м =у,, = ‚п = . 
ря и, № , > о б 
+1 
ар + Чрх, а 
| 21, 6| 0:(хьх 
ав (х)=—® ‚$ %(% к) где В=1,2,3, а=1,2,3, 
а ах. а 9, (хр, жь ) 
о" Яа [о 
п, 2, 6 
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Е @. х@,;, если В =1; 
вв (х„ хз )= 9, (х, хз Е@х@;, если В=2; 
( 


) 

з) 
х,№)е@ ха, если В=3; 

) 

) 






































9, 
(18) 
в (№,ж%), (х,ж)е @хаз, если @ =1; 
0, (вх, ) = 6, (х,ж), (м,ж)е@х@з, еслиа=2; 
9: (х,х), (х,>)е@х@, если а =3. 
Функцию @(х), соответствующую оператору О, определим так: 
з Яйх 9 
а(х)=У а +12 ба. 9 —“, ХЕО. 19 
сана] 9) 
Из соотношений (16) и (17) получаем неравенство 
оо 3 +1 о2 3 1 Чах о2 
[кр = р 9. Чу ,1| (20) 
а=| 9, “ а=| бд |2й, 6 
Заметим, что справедливо неравенство 
+Го2 02 
бу, „1 < в у 1 
В, = В. “< й 
Из этого и соотношения (18) следует: 
2 2 
в оо 3 а. 2 1 Чсх а о 
ВР" у,у |< а + | Ну | | 21 
| :) В) [= м, 6” о 
[24 а 
Будем считать, что у“ (х) — решение краевой задачи 
(вле) Ч = а =1,2,3, х, ва, (22) 
“/х, 3 2, 6 
уе =0,х, =0,х, =. 
Положим 
2 (хвьжь)= тах у” (х), х= (х,х›,ж)е @, (23) 


ХЕ, 
ге а=1,2,3. 
1 (х,ж), (№»,ж2)е@хаз, если а =1; 
в: мии |= с›(м,ж), (х,ж)е@ хе, еслиа=2; 


сз(м,%), (х,х)е@х@, если а =3. 


Используя лемму 14, приведенную в п. 4, $ 2, гл. У [8], получаем: 
к Зе 


а о? о? о? 
У ‚1 ЗС. || а, Ух ео ‚1 Г 
[2 @% 


2, 6 
Умножая обе части последнего равенства на произведение сеточных функций 9,1, 








и суммируя скалярно по 


Декартову произведению сеток @,; х@%, 











где 
а=1,2,3; 
(2,3),если @=1 
(В.д)=3 (13). ели а=2 
(1,2), ели @=3, 
получаем 
Ч х Ги] о? о? а о? 
а ол, 6] 5111 <| 9 хесь| а а? Ш (24) 
[24 


Рассмотрим теперь другую совокупность трехточечных разностных задач: 
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+1 











(али на г и =, =0, =, а =12,3; 
Ее ж= (а ео 
Пусть 
Ь, (хвьж) = тах м" (х), х=(м, 2,2 )е@, а =1,2,3, (26) 


где 
Ь (х›,%), (х,%)е хз, если а =1; 
В, (хвь хз )= В. (х,ж), (х,жз)е хе, еслиа=2; 
6: (м,^), (х,%)е@ ха, если а =3. 


Базируясь на уже упоминавшейся лемме [8], из соотношений (25) и (26) получаем оценку вида 


а! 02 02 й 02 
[45 ва [3 |} ‚ а=2,3. 
- 2, @ 


Умножим обе части последнего равенства на сеточную функцию 0., суммируя скалярно по Декартову произ- 
ведению сеток @в х 5, где 


а=1,2,3; 
(2,3), если а =1; 
(8,5) =4(1,3), если а =2; 
(1,2), если а =3, 


а# о? р 2 
[5 2}: вл у, +3 | 07) 
9 [24 


Сложим (24) и (27) и с учетом формулы (19) получим 


о 
| .|: 


получим неравенство вида 


о 


02 
| (6+ жби), а, +1 И (28) 


к 
п [Мы 


Построим функцию 
1 


т 


(29) 


Здесь 

@&=12.3, х, Е@в, Х;Е;,; 
(2,3), если @ =1; 

(8,5) =3 (1,3), если а =2; 
(1,2), если а =3. 


Тогда неравенство (28) запишется в виде 


о2 оо 
| (453) 


Следовательно, дб =1. 


Рассуждая подобным способом, придем к оценке следующего вида 


1 о? с о2 42” 
НЕ аа 30 
[= . т [ы ” — 
[24 


Неравенство (21) с учетом оценки (27) запишем 


Чех, а 


21, 6 
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оо 3 о2 о? о2 
ВО у, |< У || ву, л| + “ау +Чу |1| |. (31) 
а= 9, ы 9 “ 3 
а а 


а= тах {а(м,л.%)} (32) 


(м ‚Хо, )Е@ 








Предположим, что 


Исходя из соотношения (32), равенство (31) можно усилить и записать так: 


оо 3 о? 02 >02 
[р = у |) [ | у, +3 || ь (33) 
а=1 а я «Ха а 


Применяя лемму 17, п. 5, $ 2 гл. У [8], получаем неравенство: 


02 о2 
ыы) =. у с (34) 


(1 





где 

Та = Та (хвьх); а=1,2,3; Хх, Е @р; % Е @5; 
(2,3), если & =1; 

(8,5) =3(1,3), если а =2; 
(1,2), если а =3; 


_ 4 а (м, жж )+а, (м +1, №, 2%) | 
Я сб оо 2 


[2 (%1,%>, №) +45 (м,х гы 


2 





> 





4 
= 1 
72 РЕ пы (2,0, ) „аи (м, Е 





4 аз (м, №) аз (х,х,ж+ 
реник бильь) тах |= - 3) эм — ы 


0х. < -—№ 2 


Опираясь на неравенство (31), имеем цепочку неравенств: 


тр =? ий 
| зы =. у -и. 38 у ) &=1,2,3. (35) 
[0 © Го) 


а (2 


Будем задавать постоянные величины в соответствии с равенствами: 


к, = те. а=1,2,3. (36) 


7а +34 


Базируясь на соотношениях (35), (36), приходим к цепочке неравенств: 


02 о2 02 
аа, Ух, 1 =Аа а» Ух. 1 + (1-№л ) а Ух. 1 < 


2 





2 _# 
2 3 
кал а Па) т а = (37) 
02 — 02 
ыы зы + „| у ‚ а=12,3. 
[2 ©, 


Рассуждая аналогичным образом и учитывая неравенство (34), приходим к соотношению 


а_ 0? | 52 аз? 
| < м о +-у |1|, @=Ъ2,3. 38 
[& >. | "”% 3 | 
а [24 


Воспользуемся для преобразования неравенства (38) соотношением (33): 
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2 Е. 2 ый 
3 К [6] а о с [6] а о 
| аи ух + || + “ау +-у |1! |. 
| 9 | ый а бод е 3 а 
Используя соотношение (29), из последнего неравенства получаем оценку вида: 


оо 3 02 —о2 
[др у, У У | же. Е Ух у |: (39) 
а=1 ря “5 
[24 


<; 
Чтобы получить в неравенстве (39) наиболее точную оценку для скалярного произведения 
1 оо 
КР ЮУ, У |, 
необходимо найти минимум выражения 
[24 


(в, ников +“ | 
д 


который будет достигнут, если 
1 
(40) 





Ь р 2 
в х5) :&“=1,2,3; х Е @в; Хб Е @;; 
Ко (хв›х5) | 

(2,3), если а = 5; 


(8,5) =4(13), если а =2; 
(1,2), если а =3. 


Яа (хвьжь) э 


Неравенство (39) примет вид 


2—2 
> оо 3 С о о 
[р 1, < У | ев =] У +9, |: : 
а= Я. «3 
[24 


а 


—1 


1 
2 
+ [= : (41) 


1а 


Соответствующее выражение для определения функции 4(х) запишется следующим образом: 


1 
2|ах, а 
2, 6 











З ан Ь 
а(х)=У|-— + —% 
(*) а=1 И [= 


Из соотношения (37) и неравенства (36) получаем оценку для параметра Д: 
1 


2 

И ЕЕ 6 (хр-хь Ка (хр. р м (хр. х ] ь 

а=1,2,3, х, Е@,, х;еф,, (42) 
(2,3), если а =1 


(8,9) =3 (1,3), если а =2 
[62 если а =3. 


Ранее было показано, что имеет место оценка 0 =1. Таким образом, определить параметр @о можно с помо- 


шью следующего равенства: 
2 
@) == —= . 
УА 
Поскольку сеточный оператор А, входящий в операторную постановку задачи, задаваемую равенствами (8), 


(9), является самосопряженным, то в итерационной формуле (12) итерационные параметры ти+1, п = 0,1,...по(&) — 1 


следует для ускорения сходимости выбирать совпадающими с чебышевским набором. Тогда для числа итераций по(Е), 


(43) 


Сухинов А.И. и др Усовершенствованный попеременно-треугольный метод численного решения 





требуемого для достижения заданной ТОЧНОСТИ 85, справедлива оценка [8]: 
2 
Ала (2) 
[2 
п (Е) =. 
(=) р 


(44) 


Если № = тах { №, №, №}, 


где максимум берется по числу узлов сетки по всем трем координатным направлениям, то нетрудно показать, что 
оценка (44) принимает вид: 


и =0 муш{2 (45) 
$ 


4. Результаты численных экспериментов. Для того чтобы построить модельную задачу, будем ориентироваться на 
среднюю проницаемость нефтеносных пластов, порядок которой составляет 10" м”. Это значение взято в качестве 
базовой проницаемости. Далее базовая проницаемость умножается на неоднородный коэффициент проницаемости, 
показанный на рис. 1 (сгенерирован случайно из диапазона 0,01—100 для области 150х150 м). Таким образом, разброс 
значений проницаемости составляет 4 порядка, что соответствует наиболее сложным с вычислительной точки зрения 
постановкам задач [9]-[12], на обусловленность которых влияют как количество используемых узлов сетки, так и раз- 
брос коэффициентов сеточных уравнений. 


150 
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75 
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Рис. 1. Распределение проницаемостей пласта 

Будем считать, что: 
— толщина пласта — 10 м; 
— горизонтальные размеры области — 150х150 м; 
— граничные условия — непроницаемые стенки. 

Зададим суточные дебиты скважин, в одну из которых (нагнетательную) вода закачивается с дебитом 
300 м’/сутки, а из другой (эксплуатационной) — выкачивается нефтеводоносная смесь с отрицательным дебитом 300 
м’/сутки. Проницаемость пласта в зависимости от вертикальной координаты не меняется. Исходные данные для 
наглядности представлены в табл. 1. 





























Таблица 1 
Исходные данные к модельной задаче 
Наименование параметра Значение Единицы измерения 
Длина области по направлению ОХ 150 м 
Длина области по направлению ОУ 150 м 
Мощность пласта 10 м 
Пористость 0,2 
Проницаемость В соответствии с рис. 1 
Начальная водонасыщенность 0,1 
Дебит на нагнетательной скважине 300 м’/сут 
Дебит на эксплуатационной скважине —300 м’/сут 
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Использовалась пространственно-трехмерная сетка с общим количеством шагов: 300х300х20, шаг по времени 
составлял 0,003-0,005 сут. Сеточное уравнение для функции давления решалось двумя методами: Зейделя (МЗ) и ме- 
тодом верхней релаксации с шахматной упорядоченностью узлов (МВРШУ) [1], а также построенным в данной рабо- 
те усовершенствованным МПТМ. 

На рис. 2 представлено пространственное распределение функции водонасыщенности в нефтеводоносном 
пласте в зависимости от горизонтальных координат по истечении 20 и 50 суток, а на рис. 3 — распределение давле- 
ний. Следует отметить, что полученные результаты удовлетворяют закону материального баланса на дискретном 


уровне и качественно согласуются с результатами моделирования подобной задачи другими авторами [9]. 
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а) 6) 
Рис. 2. Распределение функции водонасыщенности: через 20 суток (а); через 50 суток (6) 
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Рис. 3. Распределение функции давления через 50 суток 
(изменяется в соответствии с палитрой от 100 атмосфер до 250 атмосфер) 


Результаты численных экспериментов демонстрируют значительное преимущество построенного варианта 
МИТМ по сравнению с широко используемыми методами — число итераций сокращается в 10—70 раз (табл. 2). 
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Таблица 2 


Количество итераций, требуемое для перехода на новый слой 























Номер временного Метод решения сеточных уравнений для функции давления 
слоя МЗ МВРШУ МПТМ 
1 37535 7232 752 
2 23347 3243 314 
3 14138 2172 196 
4 13190 1918 188 
5 12932 1863 184 




















Заключение. Для ЗО-модели фильтрации двухфазной жидкости выполнена ее дискретизация на основе мето- 
да, явного по насыщенности и неявного по давлению. Для полученной системы сеточных ЗО-уравнений построен ва- 
риант модифицированного попеременно-треугольного метода (МПТМ). Он учитывает наличие 6-образных источни- 
ков (скважин) как в структуре оператора-предобусловливателя, так и в оценке энергетической постоянной, определя- 
ющей минимально необходимое число итераций при использовании чебышевского набора итерационных параметров. 
Данный вариант МПТМ применен для решения модельной пространственно-трехмерной задачи в пласте с существен- 
но (на 4 порядка) изменяющейся проницаемостью на сетках, включающих до 200 000 узлов. В сравнении с другими 
методами данный вариант МПТМ требует меньшего числа итераций (в 10-70 раз). Таким образом, его можно рас- 
сматривать как базовый итерационный метод в алгоритме типа /МРЕ5$ — неявном по давлению и явном по насыщен- 
ности. 
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